
# =============================================================
# File: App_I_exploratory.R
# Purpose: Produces results from Appendix Section I: Exploratory Sectoral Analysis
# Paper: Foreignness as an Asset: European Carbon Regulation and the Relocation Threat among Multinational Firms
# Author: Patrick Bayer, patrick.bayer@strath.ac.uk
# Date: 5 October 2022
#
# Data: ./data.csv
#
# Technical disclaimer:
# All analyses in R version 4.1.2 (2021-11-01)
# RStudio 2022.07.1 Build 554 ("Spotted Wakerobin" Release (7872775e, 2022-07-22) for Windows)
# Windows 10 Enterprise, 64-bit
# 12th Gen Intel(R) Core(TM) i7-1270P 2.20 GHz with 32GB RAM
# =============================================================

library(lmtest)
library(sandwich)
library(MASS)
library(MatchIt)
library(cem)
library(cobalt)
library(cowplot)
library(ggpubr)
library(marginaleffects)
library(AER)
library(tidyverse)

# Load data
df <- read_csv(file = "./data.csv")

# Trim sample down to true within firm sample
# (1) Limit sample to European MNCs: no allocation data for domestically-owned operations for non-EU MNCs
# (2) Exclude MNCs with only foreign-owned operations as there is no within MNC variation

df$sample <- ifelse(df$mnc.eu==1 & df$foreign.share<1,1,0)
df <- df[df$sample==1,]
df$foreign <- as.factor(df$foreign) # set variable as categorical

# =============================================================
# APPENDIX I:  Exploratory sectoral analysis
# =============================================================

# The main text refers to exploratory sectoral results for three more industries:
# (a) paper (NACE2 code 17);
# (b) chemicals (NACE2 code 20); and
# (c) non-metallic minerals (NACE2 code 23)

# Create NACE2 indicator
df$NACE2 <- substr(df$NACE,1,nchar(df$NACE)-3) 
df$NACE2 <- ifelse(df$NACE2=="",NA,df$NACE2)

table(df$NACE2)

# Create leakage list variable
# Coded from second official European Commission Leakage List 
# Available at: https://eur-lex.europa.eu/legal-content/EN/TXT/?uri=CELEX:32014D0746

# Leakage list adopted in 2014
df$leak.list <- ifelse((
  df$NACE=="05.10" | df$NACE=="06.10" | df$NACE=="06.20" | df$NACE=="07.10" | df$NACE=="07.29" | df$NACE=="08.91" | df$NACE=="08.93" |
    df$NACE=="08.99" | df$NACE=="10.20" | df$NACE=="10.41" | df$NACE=="10.62" | df$NACE=="10.81" | df$NACE=="10.86" | df$NACE=="11.01" |
    df$NACE=="11.02" | df$NACE=="11.04" | df$NACE=="13.10" | df$NACE=="13.20" | df$NACE=="13.91" | df$NACE=="13.92" | df$NACE=="13.93" |
    df$NACE=="13.94" | df$NACE=="13.95" | df$NACE=="13.96" | df$NACE=="13.99" | df$NACE=="14.11" | df$NACE=="14.12" | df$NACE=="14.13" |
    df$NACE=="14.14" | df$NACE=="14.19" | df$NACE=="14.20" | df$NACE=="14.31" | df$NACE=="14.39" | df$NACE=="15.11" | df$NACE=="15.12" |
    df$NACE=="15.20" | df$NACE=="16.22" | df$NACE=="16.29" | df$NACE=="17.11" | df$NACE=="17.12" | df$NACE=="17.24" | df$NACE=="19.10" |    
    df$NACE=="19.20" | df$NACE=="20.12" | df$NACE=="20.13" | df$NACE=="20.14" | df$NACE=="20.15" | df$NACE=="20.16" | df$NACE=="20.17" |
    df$NACE=="20.20" | df$NACE=="20.42" | df$NACE=="20.53" | df$NACE=="20.59" | df$NACE=="20.60" | df$NACE=="21.10" | df$NACE=="21.20" |
    df$NACE=="22.11" | df$NACE=="22.19" | df$NACE=="23.11" | df$NACE=="23.13" | df$NACE=="23.14" | df$NACE=="23.19" | df$NACE=="23.20" |
    df$NACE=="23.31" | df$NACE=="23.41" | df$NACE=="23.42" | df$NACE=="23.43" | df$NACE=="23.44" | df$NACE=="23.49" | df$NACE=="23.51" |
    df$NACE=="23.52" | df$NACE=="23.70" | df$NACE=="23.91" | df$NACE=="24.10" | df$NACE=="24.20" | df$NACE=="24.31" | df$NACE=="24.41" |
    df$NACE=="24.42" | df$NACE=="24.43" | df$NACE=="24.44" | df$NACE=="24.45" | df$NACE=="24.46" | df$NACE=="25.40" | df$NACE=="25.71" |  
    df$NACE=="25.72" | df$NACE=="25.73" | df$NACE=="25.94" | df$NACE=="25.99" | df$NACE=="26.11" | df$NACE=="26.12" | df$NACE=="26.20" |
    df$NACE=="26.30" | df$NACE=="26.40" | df$NACE=="26.51" | df$NACE=="26.52" | df$NACE=="26.60" | df$NACE=="26.70" | df$NACE=="26.80" |
    df$NACE=="27.11" | df$NACE=="27.12" | df$NACE=="27.20" | df$NACE=="27.31" | df$NACE=="27.32" | df$NACE=="27.33" | df$NACE=="27.40" |
    df$NACE=="27.51" | df$NACE=="27.52" | df$NACE=="27.90" | df$NACE=="28.11" | df$NACE=="28.12" | df$NACE=="28.13" | df$NACE=="28.14" |
    df$NACE=="28.15" | df$NACE=="28.21" | df$NACE=="28.22" | df$NACE=="28.23" | df$NACE=="28.24" | df$NACE=="28.25" | df$NACE=="28.29" |
    df$NACE=="28.30" | df$NACE=="28.41" | df$NACE=="28.49" | df$NACE=="28.91" | df$NACE=="28.92" | df$NACE=="28.93" | df$NACE=="28.94" |
    df$NACE=="28.95" | df$NACE=="28.96" | df$NACE=="28.99" | df$NACE=="29.10" | df$NACE=="29.31" | df$NACE=="30.11" | df$NACE=="30.12" |
    df$NACE=="30.30" | df$NACE=="30.91" | df$NACE=="30.92" | df$NACE=="30.99" | df$NACE=="31.09" | df$NACE=="32.11" | df$NACE=="32.12" |
    df$NACE=="32.13" | df$NACE=="32.20" | df$NACE=="32.30" | df$NACE=="32.40" | df$NACE=="32.50" | df$NACE=="32.91" | df$NACE=="32.99" |
    df$NACE=="11.06" | df$NACE=="13.30" | df$NACE=="23.32" | df$NACE=="23.62" | df$NACE=="24.51" | df$NACE=="24.53" ),1,0)

# =============================================================
# Model (1): Paper industry
# =============================================================

# Share of sub-sectors on leakage list
table(df$NACE2=="17", df$leak.list)

m1 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(iso2)+as.factor(BVD), data=df[df$NACE2=="17",])
summary(m1)

# Observations
nobs(m1)
length(m1$xlevels$`as.factor(iso2)`)
length(m1$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m1)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m1, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs


# =============================================================
# Model (2): Paper industry
# =============================================================

m2 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(iso2)+as.factor(ETSsector)+as.factor(BVD), data=df[df$NACE2=="17",])

# Observations
nobs(m2)
length(m2$xlevels$`as.factor(ETSsector)`)
length(m2$xlevels$`as.factor(iso2)`)
length(m2$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m2)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m2, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs


# =============================================================
# Model (3): Paper industry
# =============================================================

# Exact matching on 4-digit NACE code and firm
df$NACE <- ifelse(df$NACE=="",NA,df$NACE) # Recode empty strings as NA 

m <- matchit(foreign ~ NACE + BVD, data=df[is.na(df$NACE)==FALSE & df$NACE2=="17",], method="exact")
df.match <- match.data(m)

# Re-run main model on matched sample
m3 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(iso2)+as.factor(ETSsector)+as.factor(BVD), data=df.match[df.match$NACE2=="17",], weights=weights)
summary(m3)

# Observations
nobs(m3)
length(m3$xlevels$`as.factor(ETSsector)`)
length(m3$xlevels$`as.factor(iso2)`)
length(m3$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m3)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m3, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs




# =============================================================
# Model (1): Chemicals industry
# =============================================================

# Share of sub-sectors on leakage list
table(df$NACE2=="20", df$leak.list)


m1 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(iso2)+as.factor(BVD), data=df[df$NACE2=="20",])
summary(m1)

# Observations
nobs(m1)
length(m1$xlevels$`as.factor(iso2)`)
length(m1$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m1)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m1, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs


# =============================================================
# Model (2): Chemicals industry
# =============================================================

m2 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(iso2)+as.factor(ETSsector)+as.factor(BVD), data=df[df$NACE2=="20",])

# Observations
nobs(m2)
length(m2$xlevels$`as.factor(ETSsector)`)
length(m2$xlevels$`as.factor(iso2)`)
length(m2$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m2)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m2, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs


# =============================================================
# Model (3): Chemicals industry
# =============================================================

# Exact matching on 4-digit NACE code and firm
df$NACE <- ifelse(df$NACE=="",NA,df$NACE) # Recode empty strings as NA 

m <- matchit(foreign ~ NACE + BVD, data=df[is.na(df$NACE)==FALSE & df$NACE2=="20",], method="exact")
df.match <- match.data(m)

# Re-run main model on matched sample
m3 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(iso2)+as.factor(ETSsector)+as.factor(BVD), data=df.match[df.match$NACE2=="20",], weights=weights)
summary(m3)

# Observations
nobs(m3)
length(m3$xlevels$`as.factor(ETSsector)`)
length(m3$xlevels$`as.factor(iso2)`)
length(m3$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m3)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m3, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs



# =============================================================
# Model (1): Non-metallic metals industry
# =============================================================

# Share of sub-sectors on leakage list
table(df$NACE2=="23" & df$NACE!="23.99", df$leak.list)


m1 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(iso2)+as.factor(BVD), data=df[df$NACE2=="23",])
summary(m1)

# Observations
nobs(m1)
length(m1$xlevels$`as.factor(iso2)`)
length(m1$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m1)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m1, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs


# =============================================================
# Model (2): Non-metallic metals industry
# =============================================================

m2 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(iso2)+as.factor(ETSsector)+as.factor(BVD), data=df[df$NACE2=="23",])

# Observations
nobs(m2)
length(m2$xlevels$`as.factor(ETSsector)`)
length(m2$xlevels$`as.factor(iso2)`)
length(m2$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m2)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m2, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs


# =============================================================
# Model (3): Non-metallic metals industry
# =============================================================

# Exact matching on 4-digit NACE code and firm
df$NACE <- ifelse(df$NACE=="",NA,df$NACE) # Recode empty strings as NA 

m <- matchit(foreign ~ NACE + BVD, data=df[is.na(df$NACE)==FALSE & df$NACE2=="23",], method="exact")
df.match <- match.data(m)

# Re-run main model on matched sample
m3 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(iso2)+as.factor(ETSsector)+as.factor(BVD), data=df.match[df.match$NACE2=="23",], weights=weights)
summary(m3)

# Observations
nobs(m3)
length(m3$xlevels$`as.factor(ETSsector)`)
length(m3$xlevels$`as.factor(iso2)`)
length(m3$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m3)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m3, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs

# =============================================================
#                       END OF CODE
# =============================================================